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1. Introduction 



Mathematical software, particularly curve and sur- 
face fitting routines, is a critical component of coordi- 
nate metrology systems. An important but difficult task 
is to assess the performance of such routines [1,2]. The 
National Institute of Standards and Technology has de- 
veloped a software package, the NIST Algorithm Test- 
ing System (ATS), that can assist in assessing the perfor- 
mance of geometry-fitting routines [3]. This system has 
been aiding in the development of a U. S. standard 
(American Society of Mechanical Engineers (ASME) 
B 89. 4. 10) for software performance evaluation of coor- 
dinate measuring systems [4]. The ATS is also the basis 
of NIST's Algorithm Testing and Evaluation Program 
for Coordinate Measuring Systems (ATEP-CMS) [5,6]. 

The ATS incorporates three core modules: a data 
generator [7] for defining and generating test data sets, 
a collection of reference algorithms which provides a 



performance baseline for fitting algorithms, and a com- 
parator for analyzing the results of fitting algorithms 
versus the reference algorithms. This paper concentrates 
on the development of highly accurate reference al- 
gorithms. 

The paper is organized as follows. We first introduce 
notation and certain key functions and derivatives that 
will be used throughout the paper. Section 2 describes 
fitting algorithms for linear geometries — planes and 
lines. Lagrange multipliers [8] are used in solving the 
constrained minimization problems, which are devel- 
oped into standard eigenvector problems. Section 3 deals 
with nonlinear geometry fitting. We use an uncon- 
strained optimization method that requires derivatives of 
appropriate distance functions. These required functions 
and derivatives are provided for the reader. Appendix A 
gives an outline of the unconstrained optimization 
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algorithm (Levenberg-Marquardt) that is modified to 
allow for normalization of fitting parameters within the 
routine. Appendix B gives, for all the geometries, the 
appropriate derivatives needed to create a valuable 
check for a local minimum. 

1.1 Notation and Preliminary Remarks 

Assume we are fitting a set of data points, {X/}, 
i = 1, 2,...,A^, that have been translated so that their cen- 
troid is the origin. Usually scalar quantities are repre- 
sented in plain type and matrix or vector quantities with 
boldface type. Other notation: 

X = {x,y,z) A point in 3-dimensional space. 

I-I The Euclidean (L2) nonn. E. g., \x\ = vx^ + y^ + z^- 

Xi = {xi,yi,Zi) The /th data point. 



x = {x,y,z) 



The centroid of the data, — ( ^Xi,^yi,^Zij- 



NV 



(Note: These and all other sums in this paper are 
taken from /= 1,2,...,N.) 



A = (A,B,C) Direction numbers that specify an orientation, A^O. 

a = {a,b,c) Direction cosines that specify an orientation. Note: 
\a\= 1. An orientation's direction numbers can be 
converted into direction cosines by: a =A/\A\. 

J The objective function, /is the sum of the squares of 

the distances from the data points to the geometry. 



M 



The N X 3 matrix containing the data points: 



Xi 


yi 


Zx 


X2 


yi 


Z2 


Xn 


yN 


Zn_ 



V The gradient of a scalar function. 

E.g.,VM.,..) = (...„M = g,g,S). 

For each geometry, we show the defining parameters, 
the equation for the orthogonal distance from a single 
data point to the geometry, the objective function, and 
a brief description of the steps in the calculation. 



2. Linear Geometries 

Linear geometries (lines and planes) are solved using 
Lagrange multipliers on a constrained minimization 
problem. Both cases reduce to a rather simple eigen- 
problem. 



2.1 Plane Fitting 

Defining parameters: 
X — a point on the plane. 
a — the direction cosines of the normal to the plane. 



Distance equation: 
di = d(Xi) = d(Xi,x,a) = a*(X/ 



-X) 



Objective function: 
J(x,a) = ^[a'(Xi 



■x)f 



Description: 

The centroid of the data must lie on the least-squares 
plane. This can be seen because V/ = at the least 
squares solution, yielding 2^ * fe — x) = 0. Multi- 
plying by l/N gives -^ fe - x ) + -^ (yt - y) 
-I-— 2]fc ~ ^) = 0. Distributing the summation gives 

a(x — x)-\- b(y — y) -\- c(z — z) = 0, which is to say 
6?(x,x,A) = 0, i.e., X lies on the least-squares plane. 
Since by assumption the data points have been trans- 
lated to the origin, and since the centroid of the data 
must be a point on the least squares plane, we can set 
x = 0. 

The direction of the fitted plane, a , can be found by 
solving the constrained minimization problem, namely, 
minimizing / subject to the constraint that \a\= 1. De- 
fine a function, G=\a\^ — 1, so that the problem is to 
minimize / subject to the constraint that G = 0. The 
method of Lagrange multipliers [8] tells us that the 
minimum occurs at a point where V/ = A VG , for some 
real number A . (Here, a,b, and c are treated as indepen- 
dent variables, since the constraint is accounted for in 
G. Therefore, V = (d/da, d/db, d/dc).) But VG = 2a, 
and V/=2(M^M)a, yielding the eigen-problem, 
(M'^M)a = Aa, referred to as the normal equations. 

This 3X3 eigenvector problem can be easily solved 
using well-established routines (e.g., Jacobi iterations 
[9]). However, we note that the eigenvectors of M^M 
are also the singular vectors (from the singular value 
decomposition) of M [9]. This allows us to gain numer- 
ical stability by applying the singular value decomposi- 
tion (SVD) to M without ever computing M^M, which 
is the method implemented in the ATS. 

Finally, we must determine how to select the correct 
eigenvector (i.e., singular vector) of the three produced 
by the SVD. The normal equations can be written as 
follows: 

2j Xi{a'Xi) = Xa 
2 yi{a'Xi) = Xb 
2 Zi{a'Xi) = Xc 
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Multiplying these three equations by a, b, and c, 
respectively, then summing the equations gives 

But the sum on the left is just the objective function, 
i.e., the sum of the squares of the residuals. Thus the 
value of A is the objective function, hence the correct 
eigenvector for the least-squares solution corresponds to 
the smallest eigenvalue. When using the SVD, we 
choose the singular vector corresponding to the smallest 
singular value, since under these conditions the singular 
values are the square roots of the eigenvalues [10]. 

2.2 Line Fitting 

Defining parameters: 
X — a point on the line. 
a — the direction cosines of the line. 



Distance equation: 

d(Xi) = d(Xi,x,a)- 



\a X (Xi — x)\ (which also equals 
x)]^) where X denotes the 



Vix/ — xP — [a'(Xi 
vector cross product. 



Objective function: 

J(x,a) = ^\a X (Xi -x)P 



Description: 

Since the data has been translated to the origin, and 
since the centroid of the data must be a point on the least 
squares plane (similar to the case above with the least 
squares plane), we set x = 0. 

The direction a can be found by following the same 
strategy as in the case of a plane. For line fitting, the 
normal equations are M^Ma = Xa just as in the case of 
plane fitting. Once again, the correct eigenvector must 
be chosen to minimize the sum-of-squares of the resid- 
uals. As shown with planes, we obtain A = ^(^ *•^)^ so 
/= — X-\-^\XiP, meaning that / is minimized when A 
is maximized. Thus the correct eigenvector choice is the 
one corresponding to the largest eigenvalue. As in the 
case of plane fitting, numerical stability is gained by 
finding the eigenvectors of M through the SVD, rather 
than by solving the normal equations. Since the singular 
values are the square roots of the eigenvalues [10], we 
choose the eigenvector corresponding to the largest sin- 
gular value. 



3. Nonlinear Geometries 

3.1 Utility Functions / and g 

The line and plane distance functions arise quite often 
in this paper, thus we define them here, calling them / 



and g respectively, giving necessary derivatives, which 
are used throughout the rest of this paper. We compute 
the nonlinear fits using unconstrained minimization al- 
gorithms, so we define the line and plane distance func- 
tions in terms of direction numbers rather than direction 
cosines. 

Let g (Xi,x , A ) denote the distance from the point, Xt , 
to the plane defined by the point, x, and the normal 
direction, a=A/\A\. The value of g is given by: 
gi = g(Xi,xA) = a ' {Xi - x) = a{Xi - x) + b{yi - j) + 
c{Zi - z). 

Let/(X/,x,A) denote the distance from the point, x,, 
to the line defined by the point, x, and the direction, 
a=AI\A\, The value of/is given by: / =/(X/,x,A) = 
\a X {Xi -x)\. That is. 



fr- 



: Vw' + V' + W', 



U - 

where v = 

w 



■-c(yi - 
■■a(zi - 
= b(xi 



y)- 

z)- 
-X) 



b{Zi - 
c{Xi - 



z) 

X) 

-y) 



This expression for/ is used because of its numerical 
stability. One should note that/ could also be expres sed 
(for derivative calculations) as / = ViX/ — xP — gl . 

Note: A , 5, and C are independent variables, whereas 
a, Z? , and c are not, because the constraint 
a^ -\- b^ -\- c^ = 1 causes a, b, and c to depend on each 
other. When dealing with the nonlinear geometries we 
treat/ and g as functions dependent on A , as opposed to 
fl, in order to use unconstrained minimization al- 
gorithms. Treating /and g as functions dependent on a 
would force us to restrict ourselves to using constrained 
minimization solvers. In the linear cases, we did solve 
constrained minimization problems. So when we differ- 
entiate with respect to A , for example, we treat a,b, and 
c all as funct ions of A, B, and C (e.g., a=A/ 
VA^ + 5^+ C^). This yields the following array of 
derivatives: 







da da da 


Va 




dA dB dC 


Vb 


= 


db db db 
dA dB dC 


y<L 




dc dc dc 
_dA dB dC_ 



\A\ 



\—a^ —ab —ac 
-ab \-b^ -be 
—ac —be \—c^ 



These algorithms normalize A at every step, so for 
simplicity of expressing derivatives, assume lAI = 1. A 
remains unconstrained; we just assume it happens to 
have unit magnitude. 
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The derivatives for/ and gi are then 



dx 



dy 



dz 



— = (x,- -x) - agi 



dB 



i = (yi-y) -^i 



;=(z,- z) - Cgi 
\agi - {Xi - x)]/fi 



dC 
dx 

^ = [bg.-(y^-y)yfi 

f^ = [cg,-(z^-z)]/fi 

M 

dA~ 

Ml 

dB~ 

Ml 
dC 



giiagi - (Xi- x)]/fi 

gilbgi - iyi-y)Vfi 
■ giVcgi - {zi - z)Vfi 



In the ATS, nonUnear geometries are fit in an uncon- 
strained manner using the Levenberg-Marquardt al- 
gorithm. The algorithm requires an initial guess as well 
as the first derivatives of the distance function. In prac- 
tice it converges quickly and accurately even with a wide 
range of initial guesses. Details of this algorithm are 
given in appendix A. Additionally, the code allows us to 
normalize the fitting parameters after every iteration of 
the algorithm. For each geometry we list the distance 
and objective functions, the appropriate derivatives, and 
the parameter normalization we use. 

3.3 Sphere Fitting 

Defining parameters: 

X — the center of the sphere. 
r — the radius of the sphere. 

Distance equation: 

d{Xi) = \Xi— x\ — r 

Objective function: 

J{x,r) = ^{\x.-x\-rf 

Normalization: 
(None) 



The above derivatives of/ are undefined when/ = 
(i.e., when Xi is on the line.) In this rarely needed case 
the gradient is given by: 

(yY^7\ VT^, VT^^, 



rVr^^, gW^v\ gW^^Y 



For cylinders, cones, and tori, the line associated with 
/is the geometry's axis. For cones the plane associated 
with g is the plane through the point, x , perpendicular to 
the axis. For tori, the plane associated with^ is the plane 
perpendicular to the axis that divides the torus in half. 

3.2 Choice of Algorithm 

Good optimization algorithms can readily be found to 
minimize the objective function, /. Usually such an 
algorithm will require an initial guess, along with par- 
tial derivatives, either of / itself or of the distance func- 
tion, di. Both sets of derivatives are given in this paper 
(those for / in the appendix). These should enable a 
reader to implement a least-squares algorithm even if the 
optimization algorithm used differs from the author's 
choice, which follows. 



Derivatives: 

ddi , ... , 

^— = - (x, - x)l\Xi -x\ 
ox 

^=-0,,-y)/|^,.-^| 

ddi / x/i I 

^— = - {Zi - z)l\Xi -x\ 
dz 



ddj_ 

dr 



= - 1 



3.4 Two-Dimensional Circle Fitting 

This case is simply the sphere fit (above) restricted to 
two-dimensions: 



Defining parameters: 

X — the X -coordinate of the center of the circle. 
y — the y -coordinate of the center of the circle. 
r — the radius of the circle. 

Distance equat ion: 



d{Xi,yi) = \/{Xi-xf-\-(yi-yf - r 

Objective function: 

/(x,y,r) = 2(V(x, - xf + (j, - yf - r)^ 



636 



Volume 103, Number 6, November-December 1998 

Journal of Research of the National Institute of Standards and Technology 



Normalization: 
(None) 



Derivatives: 

ddi 



^^ - (Xi - x)l{di + r) 

^=-(yi-y)/(di^r) 

ddi ^ _ . 
dr 



3.5 Three-Dimensional Circle Fitting 

Defining parameters: 

X — the center of the circle. 

A — the direction numbers of the normal to circle's 

plane. 
r — the radius of the circle. 



Distance eq uation: 

Objective function: 



Normalization: 
A ^A/\A\ (Here and elsewhere, "^" denotes as- 
signment of value. In this case, the 
value of A is replaced by the value 

A/IAI.) 



Derivatives: 

Mi 
dx 
ddi 
dy 
ddi 
dz 
ddi 
dA 
ddi 
dB 
ddi 

dc 

ddi 
dr 



= [gi(gi).^m).Vdi 

= Vgiigi)y+fi(fi)yVdi 

= Vgiigdz+fi(fi)zVdi 

= vgi{gi)A+fi(fi)Aydi 
= vgiigdB+mi)BVdi 

= Vgiigdc+fiifi)cVdi 
= -(fi- r)ldi 



Description: 

We use a multi-step process to accomplish 3D circle 
fitting: 
1. Compute the least-squares plane of the data. 



2. Rotate the data such that the least-squares plane is 
the x-y plane. 

3. Project the rotated data points onto the x-y plane. 

4. Compute the 2D circle fit in the x-y plane. 

5. Rotate back to the original orientation. 

6. Perform a full 3D minimization search over all the 
parameters. 

Some coordinate measuring system software pack- 
ages stop at step (5) and report the orientation, center, 
and radius as the least-squares circle in 3D. This ap- 
proach is valid when the projection onto the plane is 
done simply to compensate for measurement errors on 
points which would otherwise be coplanar. But this 
method does not in general produce the circle yielding 
the least sum-of-squares possible (even though it is usu- 
ally a good approximation.) In order to achieve the true 
3D least-squares fit, the circle computed at step (5) is 
used as an initial guess in the Levenberg-Marquardt 
algorithm, which optimizes over all the parameters 
simultaneously [step (6)]. 

Step (2) is carried out using the appropriate rotation 
matrix to rotate the direction, a, to the z -direction, 
namely. 



-ah 



\+c 



ab 



1+c 



1+c 



1- 



1+c 



-b 



If c < 0, fl is replaced with — a , thus rotating the 
direction to the minus z -direction (which is adequate for 
our purposes.) Step (5) is carried out using the appropri- 
ate rotation matrix to rotate the z -direction to the direc- 
tion, a . Namely, 



1- 



-ab 



1+c 1+c 



— ab 
T+^ 



1+c 



3.6 Cylinder Fitting 

Defining parameters: 

X — a point on the cylinder axis. 

A — the direction numbers of the cylinder axis. 

r — the radius of the cylinder. 
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Distance equation: 

d(x,)=fi- r 

Objective function: 

J(x,A,r) = J,{f,-rf 

Normalization: 
A ^AI\A\ 
X ^ (point on axis closest to origin) 



Derivatives: 

ddj 
dx 
ddj 

'dy'' 



ddj 

ac 

ddj 
dr 



(f<)y 



ddj 

dz " ^'^' 

ddi _ 

ddi _ / /■ X 
dB - ^'^^ 



(fi)c 

- 1 



Derivatives: 

dd, 
dx 
dd, 



(/;)xcost/f+(^0xsint/f 



^ (fi)ycosilj -{- (gi)ysimlj 
j^ = (fi)zCosilj-\-{gi\sirnlj 
^ = (/;)ACOst/f+ {gi)Asimlj 



= (fiXcosijj -{- (gi%simlj 
= (fiXcosijj -{- (gi%simlj 
= - 1 



ddi 
dB 

ddi 

dC 

ddi 

ds 

ddi r ' I t 

— =-fisimlj-\-g,cosilj 



3.8 Torus Fitting 

Defining parameters: 
X — the torus center. 

A — the direction numbers of the torus axis. 
r — the major radius. 
R — the minor radius. 



3.7 Cone Fitting 

Defining parameters: 

X — a point on the cone axis (not the apex). 

A — the direction numbers of the cone axis (pointing 

toward the apex). 
s — the orthogonal distance from the point, x, to the 

cone. 
t/f — the cone's apex semi-angle. 



Distance eq uation: 



d{Xi) = \/gf + {f,-rf-R 

Objective function: 

J{x^A^r^R) = ^^^/gf + {fi-rf-R^ 

Normalization: 
A 4-A/IAI 



Distance equation: 

d{Xi) =fiCosilj-\- giSimjj - s 

Objective function: 

J{x,A,s, (jj) = 2(/;cost/f + ^,sint/f - sf 

Normalization: 

A ^A/\A\ 

X ^ (point on axis closest to origin) 

(jj ^ t/f(mod 277) 

if t/f > 77 then [t/f ^ t/f(mod 77); A ^ — A ] 

77 
if t/f > — then t/f 4— 77 — t/f 

if 5<0 then [s ^ - s',A ^ - A] 



Derivatives: 

^=[g^(gi%^(fi-r)(fi%]/(d.^R) 

^=[g^(gi)y^(fi-r)(fi)y]/(d.^R) 
^=Vg.{gi)z + if-r){fi)Md. + R) 

^=Vg.{gi)A + ifi-r){fdAM{d. + R) 

dd, 
dB 

|^=U«fe)c + (^-r)(/;)c]/(^, + /?) 

-^=-(f^-rM-^R) 

ddi ^ _ . 
dR~ 



= [g^(g^)B^(f—r)(fdBV(d.^R) 
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4. Discussion 

The algorithms have been implemented in the ATS 
and have been used for 5 years by NIST and by members 
of the ASME B89.4. 10 Working Group. In general they 
have performed extremely well. They have successfully 
solved a number of difficult fitting problems that could 
not be solved by many commercial software packages 
used on Coordinate Measuring Systems (CMSs). (Some 
of the most difficult problems are cylinders or cones 
sampled over a small patch.) The ATS algorithms have 
an extremely broad range of convergence. Failure to 
converge has only been observed for pathological fitting 
problems (e.g., fitting a circle to collinear points). Spe- 
cial checks can detect most of these situations. 

The ATS algorithms are generally robust. For most 
fits, a good starting guess is not required to reach the 
global minimum. This is due, in part, to the careful 
choice of fitting parameters, the use of certain con- 
straints, and, for cylinders and cones, the technique of 
restarting a search after an initial solution is found. 



5. 



Appendix A. 
Algorithm 



The Levenberg-Marquardt 



The Levenberg-Marquardt algorithm [11] finds the 
vector,/?, that minimizes an objective function, J{p) = 
Xdfip), where di{p) is the distance from the point, X/, 
to the geometry defined by the vector of parameters, p . 
In the case of a cylinder, for example, /? = (x, A , r) = 



(x, y, z, A, B, C, r) and di{p) =fi — r. One derivation 
of the algorithm [12] begins by approximating / as a 
linear function of/?, /(/?): 

J(P)'-J(P) = J,(d,(po) + Vd,(po) pf 

where Vdt (/?o) is the gradient of dt (p ) evaluated at an 
initial guess, po. This approximation is valid within a 
certain trust region radius. The derivation then considers 
how to minimize /(/?). A search direction is calculated 
based on /(/?), and a search is made in that direction 
within the limits of the trust region radius for a point, 
/7new, such that /(/?new)</(/^o)- Whcn /7new is found, it 
becomes the new po for another iteration of the above 
process. 

The basis for the algorithm is the result that the solu- 
tion, p *, to each iteration can be expressed as 
p^=p(\)=- (Fo'^Fo + W^D)-' Fjdipo) where Fo is 
a matrix having Vdi(po) as its iih row, D is an appropri- 
ate weighting matrix, d(po) is the vector of residuals, 
di(po), and A ^ is a variable called the Levenberg- 
Marquardt parameter. This parameter can be considered 
the Lagrange multiplier for the constraint that each 
search be limited to the trust region radius. 

The Levenberg-Marquardt algorithm is presented in 
the figure below. The matrix FJFq -\- XD^D is named H 
and the vector Fjd{p^ is called v. The algorithm in- 
cludes a modification suggested by Nash [11], which is 
to use a weighting matrix defined so that D ^D is the 
identity matrix plus the diagonal ofFjFo. Nash's use of 



Input is an initial guess vector, p^ , containing the defining parameters. Returns p^ that minimizes J. 

Set A ^0.0001 

Repeat 

Decrement A ; Normalize p^ 

set U <- FoV ; V <- F,^d(p, ); J, ^ J,df(p,) 

Repeat 

Increment A 

Set iy<-l/ + /l(/ + diag(wii,M22 '•"'"««)) 

Solve the system Hx — —V 

Set />„,, <^Po+X; 7^„ <r- ^df (p^^ ) 

If converged, Set/?Q ^Normalized p^^^ ; return p^ 
Until 7^g^ < 7q or an iteration limit is reached 



If J. 

Until an iteration limit is reached 



new "^ ^0 ' Po ^ /'new 



Fig. 1. Levenberg-Marquardt algorithm. 

639 



Volume 103, Number 6, November-December 1998 

Journal of Research of the National Institute of Standards and Technology 



the identity in the definition of Z) forces H to be positive 
definite. (Note that D is never calculated.) Since H is 
also symmetric, the system Hx = — v can be reliably 
solved using the Cholesky decomposition [9]. We chose 
10 and 0.04 as factors with which to increment and 
decrement A respectively. Finally, note that we normal- 
ize the parameter vector, p , at every iteration. However, 
normalization of/7 never changes the value of the objec- 
tive function, /(p). This routine is the basis of the ATS 
implementation of the Levenberg-Marquardt algorithm. 



6. 



Appendix B. 
Correctness 



Gradient Test for 



One helpful check is to compute the gradient of the 
objective function at the computed minimum. The cor- 
rect least-squares solution yields V/ = 0. Also, these 
derivatives are important because several minimization 
algorithms require that these derivatives be included. 
The definitions of/ and g are kept the same as they were 
in Sec. 3. For every nonlinear geometry, assume the 
same definition for the objective function, /, as speci- 
fied in Sec. 3. For the linear geometries, the objective 
functions are here defined in terms of unconstrained 
parameters. 

Planes 

J{x,A) = Y,gf 



f=22^.te). 



= ^Y.8ii8dy 



dJ 



dy 

^ = ^^8i(8i)A 
^ = ^^8i(8i)B 
-^=^Yj8ii8i)c 

Lines 



J{x,A) = Jjf 






c,- - X I 



% = 2^M)c 

Spheres 

T^ = - 2^di(xi - x)/\Xi 

^= - l^diiji - y)/\x^ - x\ 

^=-2Zdi(Zi-z)/\x,-x\ 
oz 

2D Circles 

— = - 2^di(xi - x)l{di + r) 

3D Circles 

fy = 2Z[8ii8dy + (fi-rM)y] 
%. = 2Z[gi{gdz + ifi-r){fi\] 
^ = 2Z[gi{gi)A + (fi-r)(fi)A] 
^ = 2Z[gi{gi)B + (fi-r){fdB] 
^=2Z{gi{gi)c+(fi-r){fdc] 

Cylinders 

fy = 2Zifi-r)(fdy 
f^ = 2Zifi-r)if,\ 
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dJ_ 

dA'' 

dJ_ 

dB'' 

dJ_ 

dC 

%=-2l(f.-r) 



2Z{f,-r)ifdA 
■■22(fi-r)(fdc 



Cones 



2^di[(fi)jcCosil/ -{- (gi)j,siml/] 
2^di[(fi)yCosilj-\- {gi)ysmilj] 



dJ_ 
dx 
dJ 
dy'- 

ri T 

— = 2^di\ifi\cos\fj+ {gi\su\\fj^ 
az 

r\ T 

— =2^di\ifi)AC0S\fj+{gi)ASU\\fj^ 

dJ_ 

dB 

— = 2^di\{fi)cCos\fj+ {gi)csm\fj^ 

— = 2^di{- fism\fj+ gicosifj^ 
dJ 



ds 

Tori 



= 2^di[(fi)BC0sil/ -{- (gi)Bsiml/] 



= - 2Zd. 



dJ 



d, 






dz 

dJ 



|^= 2E^[^«te)c+ (/;• - r)(fdc] 
§=-25:^, 
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